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Abstract 

We present a dynamical theory of complex cooperative projects such as large en- 
gineering design or software development efforts, comprised of concurrent and interre- 
t~| , lated tasks. The model accounts for temporal fluctuations both in task performance 



effect is in agreement with empirical observations, and can be mitigated by arranging 



and in the interactions between related tasks. We show that as the system size in- 
creases, so does the average completion time. Also, for fixed system size, the dynamics 

Cn ■ of individual project realizations can exhibit large deviations from the average when 
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O^ ' fluctuations increase past a threshold, causing long delays in completion times. This 

^O ' projects in a hierarchical or modular structure. 
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1 Introduction 

One of the main challenges in a large design project or, more generally, any problem solving 
process is the coordination of the different tasks comprising that project or process |15[I5(JI 
I20| . Coordination is particularly difficult in situations where the complexity of the project 
leads to its division into concurrent, interdependent tasks whose results must be dynamically 
integrated into an overall satisfactory solution. Examples are large-scale software design 
projects, engineering design projects, and industrial research and development efforts, ft 
has long been recognized that because of the coupled nature of the component tasks such 
problem solving processes are inherently iterative in their execution j28l I5HJ . Information 
from the partial solution to a given task can trigger a chain of revisions as solutions to other, 
related tasks are modified El EH El . 

A number of formal models of cooperative processes have been proposed that predict 
the number of iterations required for completion, while suggesting optimal concurrency and 
iteration schemes which minimize the overall time or cost (e.g., [Tl ll^L 1^1251 1121 11^131)1 155] '). 
One example of these approaches is the work transformation matrix (WTM) model |43| . 
which is an extension of the design structure matrix model (see [5] for a review) . It provides 
a simple mathematical representation of design iterations that estimates the number of 
iterations and work required for a given arrangement of coupled tasks. 

Helpful as these models are, they do not account for the frequent situation in which 
projects enter a vicious cycle of continuing revisions |ll)[ I2(j| H5) resulting in budget over- 
runs, missed business opportunities and delays which, at times, can make the final solution 
obsolete 0ES1E2I- A major factor leading to these undesirable outcomes stems from the 
unpredictable fluctuations in the value that a particular unit of work on a single task brings 
to the overall project. In order to take into account these fluctuations, recent models have 
focused on individual sources of variability, including asynchronicity and random timing of 
task updates and information exchange [341 I22j , information withholding |55j , exogenous 
changes 021 l^i volatility of resource allocation (lU], behavioral choice J31, uncertainty 
of performance evaluation j^l, and the complicated landscape of performance maxima and 
minima |34|. However, since these models are not stochastic in nature, they do not account 
for the random and unpredictable nature of project dynamics, which is an intrinsic and 
unavoidable element of any complex cooperative project. 

In this paper, we present a mathematical model of iterative problem solving processes 
which explicitly incorporates the fluctuating nature of task performance and interdepen- 
dence. The model uses the formulation of the work transformation matrix to describe the 
dynamics of the project in terms of the distance to solution of the component tasks. The 
crucial fluctuating component is included in each updating step. Progress towards an over- 



all solution is thus represented by a stochastic dynamic process which captures the erratic 
nature of the evolution of real-world projects and their component tasks. We first show 
that the time to solution increases on average as the number of interactions in the project 
increase, in agreement with empirical results [TJEIEniEll ^tnd previous simulations |34[l55j . 
We also show that as the average strength of the fluctuations increases, the time to com- 
pletion increases, also in agreement with empirical studies |2'4] . We demonstrate that in a 
large project comprising many interactive tasks, a hierarchical or modular structure can, on 
average, alleviate the problem of large convergence time, as previously proposed (e.g. |I3) 
and in agreement with empirical results |^ ^J and previous studies |321 ^D d ■ 

Moreover, because of the temporal variability of the fluctuations, evolution towards solu- 
tion of any particular instance of a project can differ greatly from the average behavior. As 
a result, projects which would converge smoothly to solution in the absence of fluctuations 
can deviate significantly from this path. When the temporal variability is low, convergence 
to a solution is smooth and the finishing times are normally distributed close to the average 
value. But above a given threshold in variability, the distribution of finishing times under- 
goes a transition to a heavy-tailed log-normal form. This distribution is in agreement with 
earlier empirical results on cooperative problem solving |S| and implies possible finishing 
times far greater than the average. 

Finally, we show that the effect of fluctuations is more severe in projects which converge 
slowly to solution. This implies that hierarchical organization not only can decrease the 
average time to solution, but can also mitigate the possibly unavoidable effects of fluctuating 
interactions. 

The paper is organized as follows. In f|21we present the dynamical model and generalize 
it so as to take into account both synchronous and asynchronous updates of the problem 
solving process. In U] we study the model's average behavior and the relation between 
convergence time and project size. In f^l we consider the effects of fluctuating efforts and 
establish regimes that lead to large deviations away from the smooth convergence to the 
goal. Section |S1 summarizes our results. 

2 A Model of Group Problem Solving 

Consider a group of individuals, or even computer programs, working cooperatively towards 
the solution of a problem. An example of this process in the design space is provided by a 
large software development effort in which individuals or teams work iteratively on pieces 
of the software and pass it along to others so that whole modules and eventually the whole 
application may be completed. Alternatively, one can conceive of teams of engineers involved 



in the design of a complicated mechanism that requires tight integration of all the parts in 
order to achieve the desired goal. 

The progress of a cooperative project comprising n tasks towards a solution can be 
represented mathematically by a state vector x with n components, as in the WTM model of 
Smith and Eppinger 43 . The ith element of x represents how far task i is from completion. 
The dynamics of the system are specified by a time-dependent interaction matrix A4 which 
encapsulates the interdependencies among the tasks. The state of the overall project at time 
i + 1 is determined by 

xt+i = AjXt (1) 

where x^ is the state at time t. The initial state Xi may be taken to be a vector of I's, 
defining the scale by which we measure how much work is left on the tasks. Ideally, as the 
project proceeds, the entries of the state vector will become smaller and smaller as the tasks 
near completion. 

The interpretation of the elements of the interaction matrix is as follows (see also the 
example below). The off-diagonal elements measure how much one task's partial solution 
helps or hinders the completion of the other tasks. A positive entry signifies that one unit 
of work on task j causes (^t)ij units of rework on task i, while a zero entry means that 
tasks i and j have no direct effect on each other. In the WTM model, only positive or zero 
entries were allowed. Our theory allows for negative entries when efforts on task j's hasten 
the completion of task i. The diagonal elements of the interaction matrix Aj account for 
different rates of progress on different tasks. This is in contrast to the WTM model where 
all tasks are performed at the same rate. 

As an example, consider the following hypothetical development of an internet-based 

software client. The process has been divided into three modules: user interface, database 

front end, and network layering. It is thus modelled by a three-dimensional system as 

follows: 

work remaining on user interface 

Xt 



\ 



work remaining on database front end 
[ work remaining on network layering , 
The initial state Xq of the system is a vector of I's and the time step is one week. After a 
week, the state of the system is given by 

/ 0.65 0.2 0.05 \ f l\ / 0.9 \ 



X2 = AiXi = 



V 



0.3 0.6 0.07 
0.09 0.05 0.7 



yiy 



0.97 
0.84 



(2) 



y 0.84 J 



where Ai is a hypothetical interaction matrix which describes the dynamics over the first 
week of the project. 



Note that in this example the interdependencies between tasks slow down the overall 
progress. The tasks proceed rather quickly when taken alone (diagonal elements), but due 
to the strong interdependencies (off-diagonal elements), not much progress is made overall 
in the first week. The off-diagonal elements in this example were chosen to reflect the 
difficulties in making components of a large software effort compatible. 

Project architecture 

The matrices A4 entirely define the dynamics of the system by representing the interac- 
tions between tasks (off-diagonal elements) and the task work rates (diagonal elements). To 
a large degree, properties of the project's architecture define the elements of the interac- 
tion matrices. Such properties include the inherent interdependence of related tasks, the 
structure of the organization (hierarchical, flat, etc.), and the relationships between individ- 
uals working on tasks and communities of practice |28j . among others. Unless the project 
undergoes a major reorganization^, these interactions remain constant over the project's 
lifetime. 

Interdependencies of this type give a constant, underlying structure to the system's dy- 
namics and form a basic interaction matrix Aq . This matrix is simply the work transforma- 
tion matrix 43 , which we will refer to as the unperturbed interaction matrix or unperturbed 
WTM. The actual interaction matrices Aj in our model are created by introducing fluctua- 
tions into the elements of the unperturbed WTM. In practice, an unperturbed WTM can be 
created for a particular project by assigning numerical values to the project's design struc- 
ture matrix (DSM) . The DSM is an important tool (see 5 for a review) which encapsulates 
a project's interdependencies, usually by means of a survey of the individuals involved. For 
example, McDaniel [S2| created a WTM for the appearance design of a vehicle's interior by 
assigning each qualitative DSM entry a numerical value and averaging the values obtained 
from many sources. This process is shown in tables ^ and |2 The numerical values chosen 
for the interdependencies are of crucial importance to the system's evolution, as explained 
in sections Ol and 0] 

Fluctuations 

Interactions in a cooperative problem solving process do not remain constant on short 
time scales but vary due to a number of features, including asynchronicity of information 
exchange, exogenous changes, behavioral choice, uncertainty of performance evaluation, the 

^In the case of project reorganizations, we would introduce a second unperturbed WTM to account for 
the reorganization. The first part of the project would be modelled using the first Aq plus fluctuations, and 
the second part using the second Aq plus fluctuations. This reorganization may thus be easily accounted 
for and will not produce any interesting new dynamics beyond what is discussed in this paper. 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


1. Carpet 




W 





W 


W 











W 





2. Center console 


W 




W 








s 


w 





M 





3. Door trim panel 





W 




W 





M 








M 


w 


4. Garnish trim 


w 





M 







W 


w 





W 





5. Overhead system 


w 




























6. Instrument panel 





s 


s 


w 







w 








M 


7. Luggage trim 











w 





W 




w 


W 





8. Package tray 











w 








w 




M 





9. Seats 


w 


M 


M 


w 





w 


w 


M 




M 


10. Steering wheel 

















s 








M 





Table 1: Sample DSM characterizing interdependencies in appearance design of vehicle 
interior. Here S means strong, M medium, W weak and none. 
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Table 2: WTM for vehicle interior appearance design process [H2] created from averag- 
ing many DSM's like the one shown in tabled Here S interdependencies were given the 
value 0.15, M 0.1, and W 0.05. The diagonal elements were obtained by estimating the 
autonomous completion time for each component. 

varying relevance of the partial solutions of related tasks to each other, and other factors 
as discussed in the introduction. These factors conspire to produce fluctuations in the 
elements of the interaction matrix. The exact nature of the fluctuations in any project is 
impossible to determine in advance because of the complexity of the process and of these 
interacting factors. The problem solving effort is thus effectively a stochastic process, where 
the fluctuations introduce a random component into the dynamics. 

To see how these stochastic fluctuations affect the dynamics, consider the software de- 
velopment example introduced above. The evolution of the project in the first week was 
determined as shown in equation |5J|. During the second week, a different story unfolds, as 
the fluctuations alter the entries of the interaction matrix : 

/ 0.8 0.27 -0.05 \ / 0.9 \ / 0.9419 \ 
0.23 0.55 0.04 0.97 = 0.7725 

-0.03 0.02 0.7 



X3 = A2X2 



V 



/ 



V "•« J 



0.5524 



Note in particular that the interaction matrix has changed, so that the state vector 
changes in a different way from its evolution during the first week. During the second 
week, progress was slow on the user interface component (interaction matrix element (1,1)). 
Moreover, hindering interactions, perhaps due to a lack of communication with the front 
end developer (elements (1,2) and (2,1)), resulted in lost ground on this task. The front end 
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itself was hurt by the bad communication but programmer 2 worked hard (element (2,2)) 
and made progress. The network protocol task also made good progress as its results were 
fortuitously compatible with those of task 1 (elements (1,3), (3,1)). 

Fluctuations are reflected in our model by having the elements of the interaction matrix 
vary randomly in time with a given distribution. For simplicity, we assume that the dis- 
tribution does not vary over the lifetime of the project. It is straightforward to relax this 
assumption if the project undergoes a major reorganization as explained above. 

To focus on the effect of the fluctuations, we introduce new notation for the interaction 
matrices Aj. Recall that the off-diagonal elements of Af indicate how much one unit of 
work on a certain task helps or hinders a different task at time step t, while the diagonal 
elements represent the work rate on individual tasks. Consider a particular element of Aj. 
One contribution to this element comes from the unperturbed work transformation matrix 
(recall the abbreviation WTM), which represents some a priori prediction of the value of 
the interaction. The other contribution is from the unforeseen fluctuations. This split is 
summarized by: 

each element of unperturbed WTM , fluctuation part 

interaction matrix (constant matrix) (matrix which varies in time) 

The fluctuating part of the interaction matrix is itself a matrix of random variables that 
vary in time with certain means and variances. For each element, the fluctuation mean 
indicates how far the value of the interaction is on average from the unperturbed WTM 
value, and the variance indicates how much temporal deviation there is from this mean. 
Since the fluctuation means are constant throughout the project's lifetime, we can further 
split the fluctuation part as follows: 

each element of unperturbed WTM . fluctuation mean . 

interaction matrix (constant matrix) (constant matrix) 

. varying part of fluctuation 

(matrix of mean-zero random variables) 

The mathematical representation of this relation is 

At = Ao + (M + St), (3) 

where Aq is the unperturbed work transformation matrix, M is the matrix of mean values 
of the fluctuations, and St is a time-dependent matrix of mean-zero random variables rep- 
resenting the temporal variation of the fluctuations. Using this notation the dynamics of 
the overall problem solving process are expressed by the stochastic equation 

Xt+i = (Ao+M + St)xt. (4) 

The state vector x is a vector of random variables, and the dynamics of the problem 
solving process are characterized by the system's average value (x^) and its moments (|xi|P). 
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The monients provide a measure of how hkely the system is to be found far from its average 
behavior. 

Stopping criterion 

A final important feature of our model is the stopping criterion and the interpretation of 
the zero vector. We consider the zero vector to represent an "optimal" solution, and the 
values of the state vector as an abstract measure of the amount of work left to be done 
before a task's solution is optimal. In the evolution described by equation (^, even when 
At is constant and the dynamics are convergent, the zero vector is never reached, as the 
state vector asymptotically approaches zero as shown in figure ^ In a given project, the 
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Figure 1: Simulation of project dynamics showing the evolution of the first component of 
the state vector. This plot was created using only the unperturbed interaction matrix Aq 
shown in table [21 (no fluctuations). At left, evolution over 100 weeks; at right, a plot of 
weeks 60-100 showing the difference in completion time for two different stopping criteria 
(dotted hues). 



goal is thus reached when managers decide that the tasks' values are close enough to zero. 
This is in contrast to the WTM formulation in which the zero vector must be reached for 
work to stop. 

Our model thus implicitly contains a stopping criterion that determines when the work 
is finished. For example, the stopping criterion may be that all the elements of x are less 
than 0.05. Other, task-specific criteria are of course possible. As shown in the right panel of 
figure^ a slight change in the stopping criterion can mean a difference of weeks or months of 
in the completion time. This interpretation is quite reasonable when one considers that real- 
world projects are always terminated by managers when the solution is deemed satisfactory 
by some subjective indicators. 



Continuous time model 

The model proposed above is discrete in time, implying that all the task states are up- 
dated at once. This is applicable to situations where individuals update their information 
synchronously, such as at a weekly meeting. Conversely, in situations where information 
is continuously passed between various individuals or subtasks in asynchronous fashion, a 
continuous- time model is needed j22|. In this continuum limit one obtains the stochastic 
differential equation [2Z| 

dx= [(Ao+M-I)<it + dS(t)]x(t). (5) 

where dS = S{t)dt is a matrix of independent mean-zero Wiener processes with standard 
deviation proportional to \/dt. 

The introduction of fluctuations into the interaction matrix means that our model re- 
quires few assumptions about the nature of the project. The one assumption that is neces- 
sary, just as in the original WTM treatment, is that the project's evolution can be described 
using a linear model. This is quite reasonable because in the early stages of a project progress 
on a given task is of a general nature and can be strongly interrelated to other tasks. As 
the project progresses and nears completion, partial solutions no longer provide much new 
information and the effects of a given task on others diminishes. In other words, interdepen- 
dencies are stronger at earlier stages of the problem than at later ones. The mathematical 
expression of this weakening of the interdependencies is that the interactions are propor- 
tional to the amount of work remaining. This is a generalization of the WTM assumption 
that rework caused by task correlation is proportional to the work done on the previous 
stage. 

3 Average dynamics and project size 

Equations Q) and 101 describe the evolution of the cooperative solving process towards a 
solution. As discussed above, these equations are stochastic and describe a range of possible 
dynamics. In order to elucidate the possible behaviors, we first study the average dynamics 
of the process. As we show later, when the fluctuations in interaction strength are small, 
the evolution of the problem solving process is very likely to be close to the average. For 
larger noise, individual realizations of a given project may differ greatly from the average. 
Nevertheless, this quantity is important in studying the dependence on project organization, 
system size and interaction strengths. 

As shown in appendix ^ the matrices St do not affect the average dynamics, so the 



average state of the process evolves according to 



(xj+i) = (Ao + M)(xf) , discrete 

d(x) = (Aq + M — I)(x(i))(ii , continuous. 



(6) 



This is a simple linear system whose dynamics are well-understood. The convergence prop- 
erties are completely determined by the magnitude of the largest eigenvalue A of the average 
interaction matrix Aq-I-M. That is, if A < 1, the process will produce a satisfactory solution 
as the state vector x converges asymptotically to the zero vector; but if A > 1, the process 
will not yield a solution as the state vector diverges, or continues to increase. This behavior 
is demonstrated in figure |21 

Moreover, the value of A determines the rate of convergence (or divergence), with higher 
values of A implying slower convergence (or faster divergence). This property is explained 
below and is demonstrated in figure [3 This is the reason why the numerical values chosen 
to represent the qualitative entries of the DSM are so important, as demonstrated in the 
figure. 
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Figure 2: Simulation of average dynamics of processes with discrete dynamics and various 
values of A. This plot was created using the DSM interaction matrix for the appearance 
design of a vehicle interior (table ^ by simply varying the numerical values assigned for 
strong, medium, and weak dependencies. The diagonal entries reflect the autonomous task 
completion rate and are those of the WTM in table |2 Here M = 0. 

The above convergence properties are a fundamental result of linear algebra jlH] and have 
been noted in connection with the WTM model 031 ■ However, previous WTM treatments 
have not fully explored the elegant simplicity of this linear system, for which the asymptotic 
solution can be concisely expressed. To do this, we will assume that Aq -|- M has a simple, 
dominant eigenvalue^. Doing so we obtain (see appendix ^ for a derivation) the following 

^This condition is fulfilled by almost every reasonable interaction matrix. A simple eigenvalue has only 
one associated eigenvector, and thus its eigenspace is one- dimensional. By dominant we mean that A > jA^l 
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simple expression for the long-time solution to equations ©: 

(xf) = cuA* , discrete 

(7) 
(x(i)) = cue'-'*' "'^^* , continuous 

where the constant c is given by c = v • xq,"^ and u and v are the right and left eigenvectors, 
respectively, of the matrix Aq + M corresponding to the eigenvalue A. Expressions Q 
clearly show that the convergence rate is determined by A. It is important to note that 
the long time limit may have different meanings for different systems. If the second largest 
eigenvalue is very small compared to A, the above asymptotic limit may be reached very 
quickly, e.g. after as few as 2 or 3 iterations, as demonstrated in figure |21 above. 

The dependence of the largest eigenvector on matrix properties relates the average con- 
vergence rate to parameters such as system size and organizational structure. If the organi- 
zation of the problem solving process is flat, i.e., the level of interaction between all tasks is 
arbitrary, most of the entries of the interaction matrix are non-zero. In this case, the largest 
eigenvalue of Aq + M is given on average by (22] 

A~n^-HO(V^), (8) 

where n is the size of the matrix Aq + M and /i is the average value of its elements. If we 
further assume that the interactions are symmetric, that is, individual i's results have the 
same effect on individual j that individual j's have on individual i, a more precise result is 
available jl4| : 

X = nfi + a'^ / fi. (9) 

where a^ is the variance of the elements of Aq + M. The linear growth of A with the 
number of independent tasks (n) and the tendency of interactions to hinder progress (c) 
shows that the time to solution for a flat organizational structure grows quickly as the 
project complexity increases. The growth of A with n is demonstrated numerically in figure 

m 

In many situations, however, smaller teams work on parts of a problem which are then 
put together at a higher organizational level. From the point of view of our dynamical 

for all Ai ^7^ A. A dominant eigenvalue is necessarily real, and the dominant eigenvalue of an interaction 
matrix must be positive. There do of course exist matrices whose largest eigenvalue is not simple, as well as 
matrices with multiple largest eigenvalues, i.e. Ai, A2, • • ■ A^ such that |Ai| = IA2I = . . . = lA^j. However, it 
can be easily shown ( 35 , e.g., or see 1541 for an explicit discussion) that a nonnegative system with several 
largest eigenvalues splits into independent subsystems each with a simple dominant eigenvalue of equal or 
smaller magnitude. If some of the elements of Aq + M are negative, then there can be a dominant complex 
conjugate pair of eigenvalues. In this case the asymptotic behavior is not smooth, but oscillatory. However, 
the convergence properties are the same with A replaced by the complex absolute value |A|. 

^The case where v is perpendicular to xq is irrelevant because the values of Aq -I- M are impossible to 
exactly determine. In addition, we note that it is possible for elements of v to be very large; this situation 
corresponds to matrices whose eigenvalues are highly susceptible to fluctuations. While of great interest, 
such matrices are rare and usually not reasonable interaction matrices. This topic is quite complicated and 
beyond the scope of this paper; see )16l . or |54| for a discussion in relation to this particular problem. 
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Figure 3: Largest eigenvalue of matrices of different sizes. Samples of 1000 matrices for 
each n = 3 through 15 were generated by choosing the matrix elements randomly from a 
distribution with average /i = 0.15 and standard deviation 0.3. The solid line connects the 
average value of A for each sample, and the error bars are the standard deviations. The 
dotted line shows the analytic prediction (A) ^ nfi. 

model, this means that each task interacts only with the other tasks of the team or their 
higher order aggregates. The ensuing interaction matrix will then have blocks of nonzero 
elements and zero entries elsewhere. In this case, a slower growth of the largest eigenvalue 
with size is obtained j21| : 

A-^lnn, (10) 

implying that much larger groups can be stable when they are structured in highly clustered 
fashion, i.e. teams. 

4 The effect of fluctuations 

We now examine the effect of fluctuations on the dynamics of cooperative problem solving. 
As discussed above, these fluctuations are a key ingredient of any cooperative process and 
are due to many different factors. Recall that the fluctuations are included in the model 
by modifying the interaction matrix from the simple unperturbed WTM to the matrix 
Af = Aq + M + Si, where M is the matrix of the mean values of the fluctuations and S( is 
a matrix of mean-zero random variables representing their temporal variation. 

In the absence of fluctuations, the process is governed entirely by the unperturbed matrix 
Aq. The dynamics are thus equivalent to the average dynamics described in section |21 
with the interaction matrix Aq instead of Aq + M. That is, the process either converges 
smoothly to solution, if the largest eigenvalue of Aq is less than 1, or diverges smoothly 
away from solution if this eigenvalue is greater than 1. Moreover, for convergent processes, 
the magnitude of the largest eigenvalue of Aq determines the time to solution; processes 
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with larger eigenvalues take longer to converge. This behavior was previously discussed and 
demonstrated in figure |21 above. 

The fluctuations affect the dynamics of the process in two distinct ways. First, if the 
fluctuations have a non-zero mean, that is, they tend to help or hinder task completion, M 
is non-zero and the average dynamics are different from the unperturbed dynamics. Second, 
if the strength of the temporal variation §f of the fluctuations exceeds a certain threshold, 
the dynamics of individual project realizations may deviate substantially from the smooth 
evolution of the average dynamics, possibly causing long delays. These two effects are 
explored in sections |4 . 1 1 and [4.21 

4.1 Shifted average dynamics 

As mentioned above, average values of the fluctuations have the effect of shifting the average 
dynamics away from the unperturbed dynamics. 

The unperturbed system is given by xt+i = ApXt, and the convergence rate is determined 
by the largest eigenvalue Aq of Aq. When fluctuations are added, the average evolves 
according to 

(xt+i) = (Ao + M)(xt) 

and the average convergence rate is shifted from Aq to A, the largest eigenvalue of Aq -I- M. 
In this section for simplicity we will proceed using the discrete time formulation only; the 
results which follow are applicable in the continuous model as well. 

Myriad forms are possible for the average fluctuation matrix M. A reasonable form to 
consider is that where each element of M is proportional to the corresponding element of 
Aq by some factor to. That is, if the (2,3) element of Aq equals 0.4, indicating a strongly 
hindering correlation, we would have Af23 = OAm, while for a weakly hindering element 
such as {Ao)i2 — 0.03 we would take M12 = 0.03to,. This assumes that the relative strength 
of the correlations is not affected on average by the fluctuations. One can thus model a 
wide variety of behavior by varying m, even to have negative values. However, the most 
reasonable values for m are positive, since unforeseen events typically hinder rather than 
help the completion of the project on average. 

The average interaction matrix thus takes the form 

Ao+M = (1 + to)Ao 

and its largest eigenvalue is given simply by 

A= (1 + to)Ao. 

This value determines the average dynamics of the perturbed system, as described in section 
O In particular, for positive values of to., the average effect of the fluctuations is to increase 
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the time to solution. Moreover, if Aq < 1 but 

1 

m> 1, 

Ao 

the fluctuations cause the project to diverge on average even though an unperturbed treat- 
ment would predict convergence. This is a far more likely scenario the closer Aq is to 1, i.e., 
the slower the predicted convergence in the unperturbed system. 

For an example, please consider the average dynamics shown in figure [3 above. The 
curves in this plot were generated by varying the off-diagonal entries of Aq from table [2 
above, taking M = 0. Note that a similar plot would be obtained by keeping Ag constant 
and varying m. Assuming that the dashed line with A = 0.87 corresponds to ttt, = 0, the 
solid line could be generated by taking m — 0.103 and the dotted line by taking m — 0.207. 
In fact, given any fluctuation in this system with m > 0.149, the dynamics of the perturbed 
system will be divergent although the unperturbed system converges. Moreover, even small 
m can cause a huge increase in convergence time as demonstrated by the dashed and solid 
curves in the figure. 

Generally speaking, even if M does not take the special form studied above, a majority 
of positive entries in M indicate that the average time to solution will very likely be larger 
than the unperturbed prediction. It is a property of nonnegative matrices that increases to 
individual elements results in an increase of the largest eigenvalue 35 . Even if Aq or M 
has some negative entries, A will still be greater than Ao on average according to |H1 if the 
average entry of M is positive. 

4.2 Temporal variability 

The shift in largest eigenvalue brought about by M and described in the previous section 
can engender delays in the average dynamics, which are still described by either smooth 
convergence to, or divergence away from, solution. Far more interesting are the deviations 
from this smooth path brought about by the temporal fluctuations which result in a jagged, 
unpredictable evolution. As we now show, a high enough level of temporal variation produces 
a skewed distribution of finishing times and can result in long delays. 

Recall that the temporal variation of the fluctuations is represented by the matrix Sj, 
which is a matrix of random variables with mean zero. To be as general as possible, we 
assume that each entry of this matrix is independent of any other. Analogously to the 
previous section, we also assume for deflniteness that the typical size of the fluctuations 
is proportional to the respective average interaction values. For example, if the element 
{Aq + M)35 = 0.1, then the standard deviation of the (3,5) element of St for any t is equal 
to 0.1s, where s is a constant of proportionality. These statements can be concisely expressed 
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by the correlation rule 

{{SthjiSt')^'r) = s^[{Ao + M),j]''6u'du'Sjj', (11) 

where 6 is the Kronecker delta. 

Notice that the state of the system is determined by many random matrices multiplied 
together. This type of stochastic effect is known as multiplicative noise. One aspect of 
systems with multiplicative noise is that their distribution are approximately log-normal, 
with parameters proportional to time. For small fluctuations, the variance decreases with 
t and the log-normal distribution is very similar to a normal distribution, in which events 
far from the average are very unlikely. However, once the fluctuations exceed a certain 
threshold, the variance increases with t and the system's log-normal distribution develops a 
heavy tail at long times, implying that events far from the average are likely. 

The log-normal character of the system is shown in figures^ and below. In figure01 a 
histogram showing the work remaining on task 1 in the vehicle interior appearance design 
project (see tables^andl^l above) is shown at t = 30 for various fluctuation strengths. Each 
plot's distribution is approximately log-normal as shown in the figure by the fitted curves; 
note the differences in shape as the fluctuation strength increases. In figure |3 the time to 
solution is plotted for several levels of noise; the stopping criterion is that all components 
of the state vector be smaller than 0.05. This log-normal distribution of solution times is in 
agreement with previous results on cooperative problem solving [Hj. 

In order to understand the conditions under which fluctuations lead to large deviations 
from the average dynamics, we study the moments of the state vector. In general, the 
dynamics of systems whose second or third^ moments (|xtp) and (|xtp) are significantly 
larger than the powers of the average |(xi)p and |(x4)p are dominated by large deviations 
from the average. 

In our particular case, the moments may become orders of magnitude larger than powers 
of the average as time goes on. This is because the moments can diverge in time even as the 
average converges to zero. In this situation the problem solving process is highly susceptible 
to long delays brought about by large deviations from the average. An example of the 
deviations brought about by small and large fluctuations is shown in figure below. 

We can obtain a mathematical condition on the matrices involved which determines 
when the deviations described above will occur. To do this, recall expressions O for the 
evolution of the average state, which imply that powers of the average are given by 

|(xt)|P = cPe*Pi"-^ , discrete 

'^ ^' (12) 

|(xt)|P = cPe*J'(^-i) , continuous. 



Large higher-order moments also indicate deviations from the average but with smaller probability. 
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work remaining on task 1 at t = 30 



m = 0.02 
s = 0.2 




0.5 1 1.5 2 

work remaining on task 1 at t = 30 




12 3 4 

work remaining on task 1 at t = 30 



Figure 4: Histograms showing the distribution of the first component of x at time t = 30 over 
10,000 simulations of the vehicle interior appearance design project for various fluctuation 
strengths. The unperturbed WTM for this process is given in table |21 Recall the notation 
TO for the size of the fluctuation mean and s for the strength of the temporal variation. 
The superimposed curves are log-normal probability functions calculated using the sample 
mean and variance and normalized by the bin width and height. The curve parameters 
are jl = 0.316,0.51,0.603 and a = 0.123,0.486, 1.019 respectively by increasing fluctuation 
strength. 
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m = 0.01; s = 0.1 
m = 0.02; s - 0.2 
m = 0.04; s = 0.4 




50 100 150 200 250 300 350 400 

finishing time 

Figure 5: Distribution of finishing times for various fluctuation strengths over 10,000 
simulations of the vehicle interior appearance design project. The curve parameters are 
fi — 109.9, 142.61, 191.07 and a — .0831, .179, .3238 respectively for increasing fluctuation 
strength. 
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Figure 6: Example of dynamics for various fluctuation strengths. At left, the fluctuation 
causes a delay in task 1. At right, task 1 proceeds more quickly than an the unperturbed 
treatment would indicate. Recall that the parameter s from (|ll|l measures the strength of 
the fluctuations. 

for large t. Here we are using the notation A for the largest eigenvalue of Aq + M and we 
have normalized the right eigenvector u to have length 1, as is customary. In general, the 
moment evolution at large t can also be expressed in this form: 



ii^tn^c- 



— /Pp^-^p 



(13) 



(or (|x(i)|P) for the continuous case) where Lp is the pth moment Lyapunov exponent and 
is always greater than the corresponding exponent in (|12|l . The interesting case described 
above occurs when A < 1, so that the average dynamics are convergent, but Lp > 1 for 
p = 2 or 3 so that the second or third moments diverge. 

The Lyapunov exponents are generally difflcult to calculate. However, jHl] gives a mea- 
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sure of the difference between the Lyapunov exponents and the corresponding exponents in 
the powers of the average: 

e^='^Y.'^n{Ao+M\,]^ul (14) 

Large values of e^ indicate that Lp is significantly larger than pin A. In this case, even a 
process which converges rapidly on average, i.e. with a low A, may undergo large deviations 
and encounter delays. When e^ is small, deviations are unlikely unless the convergence of 
the average is very slow, i.e. A is close to 1. Since the variance of the (i,j) element of 
the fluctuations is given by s^[Ao + M)y]^, it is clear from equation H14|) that the effect of 
fluctuations is increased as their magnitude increases, as shown in figure |^ 

To show analytically how e^ affects the dynamics, we consider cases in which the Lya- 
punov exponent in equation H13|l can be approximated calculated. When the evolution is 
discrete and e^ small, the moments evolve according to [SJ 

(IxfP) = c^e*'^^'"'^+'^ ' , second moment 

. (15) 

(|xtp) = g3gt(3inA+3£ ) ^ i\md moment. 

These expressions show that the moments will diverge as e^ increases. 

In the case in which the elements of the fluctuations all have the same variance iP' , |K1] 
provides a measure of the critical strength of this variance above which deviations are likely: 

' (16) 



n - 






When the variance of the fluctuations exceeds b^, the second moment of the system will 
diverge. By applying the relation IjHJl we can find a corresponding critical value s^, and thus 
determine the stability diagram of the project dynamics. The critical values are plotted 
for various n in figure 13 When A < 1 but the fiuctuation strength is larger than s^ or 
&^, deviations from the average are likely and although the process converges on average, 
significant delays are likely to occur in particular realizations of the project. 

In the continuous case, analytic results are only available for a two-dimensional system 
(two correlated tasks) and small noise. In this case, the moments evolve according to [^ 

(|x,|P) = c^exp (i[p(A- 1) +pV/2+p7i + pO{e^) + 0{p^)]) 

where 71 is a complicated function of the noise and the eigenvectors u and v and O(e^) 
means an unspecified expression whose magnitude is on the order of e^. 

5 Conclusion 

In this paper we presented a dynamical theory of cooperative problem solving projects 
consisting of any number of interrelated tasks. The model incorporates the unavoidable 
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Figure 7: Stability diagram of a typical system with n = 5. Recall that b is the absolute 
strength of the fluctuation, while s is the strength proportional to the average size of the 
entries of Aq + M. 

fluctuating nature of the interactions in such projects, which reflect the natural variations 
in the amount of progress made on each task, as well as the varying relevance of information 
from the partial solutions of each task to any other. While simple, the theory captures 
key aspects of the dynamics, such as the diminishing value of interactions as the tasks near 
completion, as well as allowing for a satisficing decision of when the project is completed. 
It also captures the fact that projects with too many tasks do not converge to a desirable 
solution at all. 

By analyzing the stochastic dynamics we showed that as the strength of the fluctuations 
in task effectiveness increases past a certain threshold, the convergence of a typical process 
to a solution is very likely to deviate from what an average treatment would indicate. In 
order words, single problem solving instances are likely to deviate significantly from what 
one would expect if fluctuations were ignored. In particular, a cooperative problem solving 
process may be pushed far from its goal, resulting in a significant delay in the time to 
solution. While this effect is most severe in problems whose average convergence is slow, 
large fluctuations can also disrupt the smooth convergence towards a solution in any project. 
This prediction is in agreement with numerous empirical results. 

Moreover, over appropriate time intervals the distribution of performance among tasks is 
log-normally distributed, as shown in previous work on cooperative problem solving efforts. 
The appearance of a heavy tailed distribution above a threshold value in task variability 
implies that repeated observations of the evolution of a project will tend to markedly differ 
from the average dynamics that other models have analyzed. In other words while there 
might be instances where time to completion seems to be in agreement with an estimate of 
the average lifetime of the project, on other occasions project duration will be dramatically 
longer, with the concomitant aggravation of not being able to predict when such long delays 
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will take place. 

These results suggest that in order to have a predictable outcome in large cooperative 
projects, effort should be made to organize cooperative tasks in a way that minimizes both 
the instabilities implied by large unstructured groups, and the fluctuating nature of their 
contributions. This can be achieved by structuring tasks in such a modular or hierarchical 
way so that they effectively interact with very few others. Furthermore, given the stochastic 
nature of such processes, one possible mitigating design factor is to keep fluctuations in 
performance constrained to variations far below their threshold value. 

Given the importance of large teams in the design and solution of technical problems 
and the increasing need for production agility, the study of the dynamics of group problem 
solving offers new insights into the organization of teams and their efficiency at solving 
complex tasks. 
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Appendices 

A Fluctuations have no effect on the average 

As mentioned above, the average dynamics of the system are not affected by the fluctuations, 
fn order to see this, notice that in the discrete case the system evolves according to xt = 
(C + Bt)xt_i. Since the fluctuations are Markovian, i.e. the components of the matrix Bt 
are independent of the state at previous times, 

(x,) = [C + (Bi)](x,_i) 

and we are taking the expected value of Bt entry by entry, ((i?T-)ij) = for all r, i and j, 
we have (xt) = C(xt_i) and, extrapolating back to t = 0, we obtain the discrete version of 
equation 10. Similarly, in the continuous case formal integration gives 

x(i) ^ I (C - I)x(s)d.s + I dWx{s). 
Jo Jo 

Since dW = 'B(t)dt and the elements of B have mean 0, 

/ /"dW^(x,s)\ =0 

for any function 0(x, s) and the expected value is given simply by 

(x(i))= /(A-I)(x(s))ds, 
Jo 

which is equivalent to the continuous version of equation . 

B Dynamics of a simple linear system 

For any dynamical system specified by the equations 

(xf+i) = C(xf) , discrete 

(i(x) = C(x(t))(ii , continuous 

the state at any time may be expressed in terms of the initial state by writing 

(xf) = C*xo , discrete 

(x(i)) = e(*-'"')*x(0) , continuous. 

In the long time limit, the matrix product C* and the matrix exponential e'*^^''* in these 
equations are dominated by the contribution from the largest eigenvalue for large t. To 
see this, we diagonalize^ C into the form C = PAP^^, where A is the diagonal matrix of 

''It is not necessary to consider defective (non-diagonalizable) matrices in our case because it is impossible 
to exactly determine tiie values of the elements of the interaction matrices. 
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eigenvalues of C. This simplifies equations H17|) because 

C* = PA*P^ 



Ai 



A| 



where the {Ai} are the eigenvalues of C, and 

p(A2-l)t 



„(A„-l)t 



p-1. 



Since A is the largest eigenvalue, when t is large A* and e'^'^"'^-'* are much, much greater than 
all other terms in the above matrices. The above equations are thus dominated by A and its 
eigenvectors u (corresponding column of P) and v (corresponding row of P~^). Expanding 
the matrix product with all entries of A neglected except the ones with A, we obtain ((TJ. 

C Log-normal distribution 

The log-normal probability distribution describes a random variable whose logarithm is 
normally distributed. The probability density function is 

1 



P(x) 



(In a;- In m) 



(18) 



where /i and a are the mean and standard deviation of the normal distribution that describes 
In a;. The pth moment of the log-normal distribution is given by 



Tf (J 

{xP) = exp {p In /i + ^— ) • 



(19) 



In our application, a^ and In/i are proportional to t. For large a^ the moments diverge for 
large enough p, causing large deviations from the average. 
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